library(SDMTools)#load necessary libraries
work.dir = 'E:/SEABat/'; setwd(work.dir) #define & set the working directory

###01.create a vector of species names
files=list.files(pattern='asc') #get a list of files in your working directory that include 'asc'

species=NULL #create the object name of your vector
spp=NULL #create a temporary object name for elements of your vector

for(file in files) { cat(file,'\n') #loop through the files, and for each file, complete the following:
spp=strsplit(file,'\\.') #split your file name by '.'  the two resulting elements will be 'Rhinolophus_sp' and 'asc'
spp=spp[[1]][1] #select the first of your two elements and call it 'spp'.  spp will be 'Rhinolophus_sp'
if(length(species)==0){ #if species is still null, ie no species in your vector yet
	species=spp #then add the first species
	} else { #once you have one species already in the vector
	species=c(species,spp)} #bind each new species to it
}

###cycle through your vector of species, and find their threshold, create a png, and determine species richness
#setup some plot parameters
legend.pnts = cbind(c(93,91,91,93),c(2,2,7.5,7.5))

richness=NULL
for(spp in species) {cat(spp,'\n')#loop through your species vector
###02.find the threshold
line=readLines(paste(spp,'.html',sep=''))[12] #read in line 12 of your html file
line=strsplit(line,'<td>')#split line 12 by <td>
line=strsplit(line[[1]][23],'<')#split the 23rd element by <
threshold=line[[1]][1]#your threshold is the first element
threshold=as.numeric(threshold) #convert your threshold to a numeric value


###03.create an image
#setup some plot parameters
bins = seq(0,1,length=101); bins = cut(threshold,bins,labels=FALSE) # get the threshold bin for cols
cols = c(rep('gray86',bins),colorRampPalette(c("wheat","tan","forestgreen","darkgreen"))(100)[bins:100])

#readin the necessary asciis....
bat.asc = read.asc(paste(work.dir,spp,'.asc',sep=''))

#modify above plotting to create single for current
png(paste(gsub('_',' ',(strsplit(spp,'//.'))),'.png',sep=''),width=nrow(bat.asc),height=ncol(bat.asc),units='px',res=300,pointsize=6) #start the png file
	par(mar=c(0,0,0,0)+0.1) #set the margins of the plot
	image(bat.asc,axes=FALSE,ann=FALSE,zlim=c(0,1),col=cols) #image the first map and add associated text & lat/lon info & legend
	legend.gradient(legend.pnts,col=cols, c(0,1), title='Suitability',cex=2)
dev.off()

###create binary ascii and store in memory

bat.asc[which(is.finite(bat.asc) & bat.asc>threshold)] = 1 #set all parts of the species distribuiton to 1
bat.asc[which(bat.asc<1)]=0
	if (is.null(richness)) {richness = bat.asc;  } #define the output & set everything = 0
	richness = richness + bat.asc #append the species-specific distribution info

}
write.asc.gz(richness, "bat_richness.asc") #write out the richness data
base.asc=richness
base.asc[which(is.finite(base.asc))]=0
richness[which(richness==0)]=NA
max.richness = max(richness, na.rm=TRUE)
min.richness=min(richness,na.rm=TRUE)
#threshold=1
#bins = seq(0,1,length=101); bins = cut(threshold,bins,labels=FALSE) # get the threshold bin for cols
#cols = c(rep('gray86',bins),colorRampPalette(c("wheat","tan","forestgreen","darkgreen"))(100)[bins:100])

cols = colorRampPalette(c("tan","forestgreen"))(100)

 #get the max richness value
png('bat_richness.png', width=nrow(bat.asc),height=ncol(bat.asc),, units='px', res=300, pointsize=5, bg='white') #start the plot
	par(mar=c(0,0,0,0)+0.1) #remove any plot margins
	image(base.asc,ann=FALSE,axes=FALSE, col='gray86')
	image(richness, ann=FALSE,axes=FALSE,col=cols, zlim=c(min.richness,max.richness),add=TRUE) #plot the richness
	legend.gradient(legend.pnts,cols=cols,limits=c(0,max.richness), title='Richness', cex=2)
dev.off()
